Power-law distributions in a trapped ion interacting with a classical buffer gas 
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Classical collisions with an ideal gas generate non-Maxwellian distribution functions for a single 
ion in a radio frequency ion trap. The distributions have power-law tails whose exponent depends 
on the ratio of buffer gas to ion mass. This provides a statistical explanation for the previously 
observed transition from cooling to heating. Monte Carlo results approximate a Tsallis distribution 
over a wide range of parameters and have ab initio agreement with experiment. 



The behavior of ions in the collision-free regime of a ra- 
dio frequency ion trap is well understood. Laser-cooling 
and the properties of the quantum mechanical ground 
state [H have been examined in great detail. It is there- 
fore surprizing that the more accessible regime of ions 
cooled by buffer gas collisions has never been thoroughly 
analyzed. Dehmelt hrst showed in 1968 that collisions 
with neutral buffer gases could either cool or heat the 
ion Q, depending on their relative masses. His theory, 
still widely accepted today, hypothesized that the recoil 
from a collision with a heavy neutral atom heated the 
ion by disrupting its response to the rf field (i.e. micro- 
motion). The theory relied on the "pseudopotential" or 
time-averaging approximation and did not address the 
statistics of the ion's distribution function, which were 
assumed to be non-critical. Subsequent workers have in- 
troduced Maxwell-Boltzmann (MB) statistics in several 
ways, for example, by assuming a Gaussian velocity dis- 
tribution function [3| or by hypothesizing Gaussian ran- 
dom noise in a Langevin equation Numerical work 
has shown apparent instability in individual trajectories 
without addressing the statistics 0. 

In this Letter we compute the ion's distribution func- 
tion using a combination of Monte Carlo and analytic 
methods. Our results show that the distribution is not in 
general Gaussian and does not follow MB statistics. Col- 
lisions with heavy neutrals give the distribution function 
power-law tails E~" in place of the Gaussian equivalent 
exp(— E I kT) . Here E is a time-averaged "pseudoenergy" 
which is not conserved during collisions^. In previous 
work the instability was thought to arise from a positive 
heating rate dE/dt >0, leading to exponential runaway. 
Our results lead to a different picture in which station- 
ary power-law tails lead to a small but constant rate of 
ion loss. This leads to orders-of-magnitude differences in 
predicted ion lifetimes. 

In the last decade non-Gaussian statistics have been 
studied in many different contexts [fjj]. In atomic physics 
such distributions have been observed primarily in laser- 
cooled atoms operating near the quantum-mechanical 
ground state. Subrecoil laser cooling has been shown 
to obey Levy statistics with a nonstationary distri- 
bution and nonergodicity H, [§]. Levy walks, anomalous 
diffusion [IgT|. and a tuneable Tsallis distribution 11 1 have 



been observed in optical lattices. The present case is 
remarkable in that it shows non-Gaussian statistics in 
a classical system not far removed from the ideal gases 
originally studied by Maxwell and Boltzmann. 

An additional motivation for understanding collisional 
heating is that ion traps have recently come into use 
as probes of collision physics. Hybrid traps comprising 
both ion traps and MOT's (magneto-optical trap) have 
been constructed yielding results for ch arge exchange 
cross-sections 12j and radiative lifetimes [13(. The the- 
ory of ultra-cold atom-ion collisions has been developed 
SEE [l6| and novel effects of an ion in a BEC (Bose- 
Einstein condensate) have been predicted [13] • Room 
temperature buffer gases have been used to study molec- 
ular ions [IH and the properties of multipole traps [l9j ]. 
This is in addition to the more traditional use of buffer 
gase cooling in trace element detection 2(J 21 1, the trap- 
ping of radioactive ion beams [B| , and several other appli- 
cations. It is necessary to understand collisional heating 
to disentangle the effect of the trap fields from the colli- 
sion physics. 

Previous Monte Carlo work has used numerical inte- 
gration to compute the ion's motion between collisions. 
This is neither fast nor accurate enough for the total 
of sa 10 10 collisions needed to compute the distribution 
function. Instead we use the classical time-evolution ma- 
trix M of the ion 
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to propagate the position and velocity of the ion from one 
collision to the next. Consider first the case of a simple 
harmonic oscillator, which obeys x(t) + (3 2 x(t) = 0. The 
time-evolution matrix S of this system is 
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This is a special case of the general solution[2 
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where P and Q are the two linearly independent solutions 
of a second order linear differential equation, R and S 
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TABLE I: Tsallis parameters n and qr fit from Fig. 1 



Distance x from Trap Center (cm) 



FIG. 1: Monte Carlo distributions for a sine le 136 Ba+ ion 
cooled by six different buffer gases at 300K ranging from 
m_e=4 (left) to m_s=200 (right). Note the evolution from 
Gaussian to power-law (straight line) as the mass increases. 
The solid lines are Tsallis functions Eq. 7 with fixed a = 
0.0f85 cm and the exponents of Table I. 



are the respective time derivatives, and the Wronskian 
D=SiPi-QiRi. The rf ion trap obeys aMathieu equation 
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for which P and Q are given by the Fourier solutions [H 0] 
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where Qi has sin[(/3 + 2m)ti)\ in place of the cosine. Here 
q= 2eVo /mi£l 2 r 2 : a= 4e£/o / niiQ 2 r 2 , E and Uq are the rf 
and dc applied potentials, ro is the trap radius, and the 
unit of time is 2/f2, where fi is the angular frequency of 
the applied rf. The coefficients up to and including C±s 
are evaluated to 24 bit accuracy (5 x 10 -8 error) with 
a recursive routine. Error propagation has been tested 
with the identity M{t N ,t x ) = llf^ 1 M(t i+1 ,U). For N 
up to 10 6 the discrepancy < 1 x 10~ 5 for randomly chosen 
times U. 

The distribution function is computed by Monte Carlo 
averaging over six random variables for each collision: the 
time U, the center-of-mass angles Qi and tfi, and the three 
random velocities vf ,z;^,and vf of the buffer gas, which 
obey a Maxell-Boltzmann distribution at temperature T, 
where T=100K, 300K, or 1000K. We assume a linear ion 
trap with transverse rf confinement and a DC potential 
along the z-axis, represented by the Mathieu matrices 
and My, where q y = —q x , and static harmonic oscillator 
matrix S z as in Eq. 2 above. We combine the three 
matrices for the x,y, and z axes into a single 6 by 6 matrix 



Buffer Gas 


mi /m fl 


n 


qr 


He 


34.5 


>60 


1.03 


Ar 


3.40 


8.2 


1.12 


Kr 


1.70 


3.8 


1.26 


Xe 


f.O 


1.98 


1.51 


170 


0.80 


1.50 


1.80 


200 


0.68 


1.15 


1.87 
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where M has M Xl M y and S z along the diagonal, 
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where, for example, r x is the column vector (x, x) 
lisions are represented by a matrix C which leaves the 
coordinates unchanged but transforms the velocities ac- 
cording to a hard sphere (isotropic in the center-of-mass) 
collision model. The ion-neutral atom collisions are mod- 
eled by Langevin scattering in which the cross-section 
(7 oc 1/v so the time between collisions is a random vari- 
able independent of the relative velocity v. An ensem- 
ble typically consists of 10 6 trials each containing 500 
to 50,000 collisions. All distributions for >500 collisions 
agree within statistics. 

Distribution functions for a single 136 Ba + ion at q=0.1 
are shown in Fig. 1. Six different buffer gases with masses 
tub = 4, 40, 84, 136, 170 and 200 amu have been as- 
sumed, corresponding to the noble gases He, Ar, Kr, Xe 
as in a recent experiment 2^, 21 1 and to two hypothet- 
ical heavier gases. The distribution for 136 Ba + in He 
is a good fit to an MB distribution with a classical a w 
yj2kT jmiuS 2 , where ui is the secular frequency {3fl/2. All 
of the other gases show non-Gaussian distributions which 
develop broad power-law tails as the mass increases. The 
four heaviest gases fit a power-law x~ 2a with good x 2 
over at least 3 orders of magnitude. For = 84, 136, 
170, and 200 the best fits are a =3.2, 1.98, 1.5, and 1.17 
respectively. A typical fit will have x 2 / n < 1-1 f° r n > 100 
degrees of freedom. 

In the absence of an analytic theory, the data has been 
fit to a Tsallis function T(x/<j,n) 0] 



T(x/(J, n) 
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which is a generalization of the Gaussian. For mooT 
reduces to a Gaussian while for small n it has power-law 
tails of the form (x/ct)" 2 ™. The exponent n is related 
to the more familiar "entropic" Tsallis parameter qr by 
qT=l+l/n . The Tsallis function arises in the theory of 
nonextensive entropy Q but at present we treat it empir- 
ically. Table 1 shows the value of n extracted from fitting 
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the distributions to Eq. 7, where a was held constant for 
all tub and To normalizes the distribution to unity. The 
fit is qualitative since the x 2 is poor due to systematic 
deviation near the origin, where the standard deviation 
< 0.3%. Nevertheless Fig. 1 shows good agreement over 
a factor of 10 5 in probability density and a factor of 50 
in buffer gas mass. The value of the Tsallis exponent n 
is close to the value of a extracted from the power law fit 
above. Similar data at 100K and 1000K give comparable 
fits with the same a scaled by VT. 

The usefullness of the Tsallis function is that it shows 
that n and a act independently of each other, to first or- 
der. In the light gas limit — ► 0, the ion has a Gaussian 
distribution with a temperature T equal to that of the 
buffer gas, so that a ~ \JlkT /rniui 2 . As mg increases, 
a changes very slowly, so that the distribution retains 
a Gaussian-like core of constant width as the power-law 
tails get stronger. This indicates that the increase in the 
mean energy of the ion is not the cause of ion loss. A 
three-parameter fit, in which er, n, and the normalization 
are optimized for each value of mg , shows a weak depen- 
dence of a on mj/iEs. For example, the best value of 
<7=0.0f75 cm at hib=4 rises to <r=0.022 cm at hib=200, 
an increase of only 26 % for a 50-fold decrease of mj/rnB. 
Similarly, changing the temperature of the buffer gas does 
not alter the power-law exponent. To generalize further, 
the Tsallis exponent n is approximated by the simple re- 
lation n sa 2m/ /m^ , which is accurate in the exponent 
to about ±20%. 

The Monte Carlo also computes the ion lifetime. The 
ion is started at the origin with zero energy and is prop- 
agated through i collisions, until r^= ^ xf + yf > r , the 
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where 



trap radius. In general the ion lifetime r oc 
n is the Tsallis exponent of Table 1. However, since the 
trap depth U oc r§ (for constant q), it is more general to 
plot r versus U, which yields r cx U™ as shown in Fig. 2. 
Interestingly, r is not sensitive to initial conditions and 
an ion starting with an energy 1 eV has r only slightly 
shorter than with zero energy. This is because most of 
the hot ions equilibrate to 300K in a few dozen collisions. 
It is only in a very large ensemble (10 6 trials in Fig. 1 ) 
or a very large number of collisions (N=5 x 10 5 in Fig. 
2) that extreme values of r^ are reached. In contrast to 
the exponential runaway model, Fig. 2 suggests that ion 
traps may be designed to achieve a specific ion lifetime. 

Power-law tails dominate the trap stability whenever 
a stationary distribution function exists. However, when 
the Tsallis exponent falls below 1, which occurs for raj > 
1.55 to,/, the distribution becomes time-dependent, as in 
a Levy distribution ,7j, the mean values diverge, and the 
ion's energy increases with each collision. In this case 
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exponential runaway occurs as originally suggested in 0] . 

The Monte Carlo results have no free parameters and 
the predicted ion lifetimes agree with the results of a re- 
cent experiment 2^, 2l| in which a single Ba + was con- 
fined in a trap of radius Rt=0.26 cm and q=0.52. Stable 




Trap Depth (eV) 



FIG. 2: Predicted power-law ion lifetime r oc U n where U is 
the trap depth, n is the Tsallis exponent of Table 1, and the 
parameters correspond to Fig. 1 



trapping was observed for He gas, Ar gas was measured 
to give an ion lifetime of 50-100 sec, while Kr and Xe 
had lifetimes too short to measure (< 5 sec). When the 
above Rt and q are input to the Monte Carlo it yields 
lifetimes of 45 sec for Ar and < 0.1 sec for Kr and Xe. 

It remains to provide a physical explanation of these 
results. An analytic expansion of the Mathieu matrix Eq. 
3 and Eq. 6 shows that the power-law tails are a result of 
a multiplicative random process 24j, i.e., products of ran- 
dom variables. These may be contrasted with the better- 
known additive random processes (sums of random vari- 
ables) , which obey the central limit theorem and produce 
Gaussian statistics. Multiplicative random products are 
not well understood but they do not in general lead to 
Gaussian distributions. Multiplicative fluctuations have 
recently been studied in a Langevin equation (25| and 
have been shown to lead to a tuneable Tsallis distribu- 
tion. 

A model for multiplicative fluctuations can be derived 
by expanding the Mathieu matrix Eq. 3 to first order in 
q and substituting the result in Eq. 6. In the same limit 
used by Dchmclt|2] one can show[23j that the matrix 
product Eq. 6 can be approximated by a product of 
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in trace atom detection [U[2l| and in trapping radioac- 
tive ion beams It is also necessary for understanding 
recent ion trap collision experiments ref[ll-17], since the 
non-Gaussian distribution function can alter their inter- 
pretation. 



FIG. 3: A heating and cooling diagram for a single collsion in 
the case mt, = rrn or a = in Eq. 8. Heating occurs for R> 1 
and cooling for R< 1. The explanation given in the text does 
not depend on the exact form of this result and requires only 
that there be a small phase space volume /3 < 1 for R > 1. 



Here (pi and (p l m are the phases of the secular motion 
(" macromotion" ) and driven rf oscillations ("micromo- 
tion") at the time of the i-th collision and a — (mi — 
mB)/(jnj + tub) is a recoil parameter. Multiplicative 
random products tend to be dominated by rare events of 
large amplitude [24]. In the present case these events can 
be identified as N consecutive heating collisions without 
an intervening cooling collision. Consider a volume of 
phase space (3 < 1 around the heating maximum 
in Fig. 3. N consecutive collisions will give an amplitude 
of 3 N / 2 with a probability of f3 N . This provides a mech- 
anism for the power-law tails leading to ion loss. If the 
trap were of infinite size, cooling collisions, which occupy 
most of the phase space, would eventually return the ion 
to the origin. In this sense, the ion loss is due to a non- 
Gaussian fluctuation rather than to heating. In the light 
mass limit a — > 1 Gaussian statistics return since each 
term in the product is near unity, R= l + e((p s , (p m ) where 
e << 1. The product reduces to a sum 1 + ^ e(cp s , <p m ), 
the fluctuations become additive, and the central limit 
theorem applies. 

This work has several implications. For statistical me- 
chanics it provides a simple, classical system which shows 
tuneable non-Gaussian statistics. For trapping and cool- 
ing experiments it shows how traps may be engineered for 
a specific ion lifetime, as in Fig. 2. This should be useful 
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